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Abstract —A framework is introduced to integrate renewable 
energy sources (RES) and dynamic pricing capabilities of the 
smart grid into beamforming designs for coordinated multi-point 
(CoMP) downlink communication systems. To this end, novel 
models are put forth to account for harvesting, storage of nondis- 
patchable RES, time-varying energy pricing, as well as stochastic 
wireless channels. Building on these models, robust energy 
management and transmit-beamforming designs are developed 
to minimize the worst-case energy cost subject to the worst- 
case user QoS guarantees for the CoMP downlink. Leveraging 
pertinent tools, this task is formulated as a convex problem. 
A Lagrange dual based subgradient iteration is then employed 
to find the desired optimal energy-management strategy and 
transmit-beamforming vectors. Numerical results are provided 
to demonstrate the merits of the proposed robust designs. 

Index Terms —Smart grid, renewable energy, downlink beam- 
forming, CoMP systems, robust optimization. 


Nomenclature 
A. Indices, numbers, and sets 

T, t Number and index of time slots. 

/, i Number and index of distributed BSs. 

K, k Number and index of mobile users. 

S, s Number and index of sub-horizons. 
j Iteration index of the dual subgradient ascent. 

T Set of the total scheduling horizon. 

Ti^s Sub-horizon s for BS i. 

K. Set of mobile users. 

Ti\. Channel uncertainty set of user k in time slot t. 
£i Uncertainty set of the energy harvested at BS i. 

£l, £l Polyhedral and ellipsoidal cases of £i. 
i Iteration index of the Bundle method. 


B. Constants 

cr^ Variance of channel additive noise. 
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a*, /3*; Ip*, (p* 


Radius of the channel uncertainty region 
for user k in slot t. 

Target signal-to-interference-plus-noise- 
ratio (SINR) of user k. 

Lower and upper bounds of the renewable 
energy at BS i in slot t. 

Lower and upper bounds of the renewable 
energy at BS i over the sth sub-horizon. 
Initial battery energy level of BS i. 
Battery capacity of BS i. 

Minimum and maximum (dis-)charging 
power at BS i. 

Battery efficiency at BS i. 

Fixed energy consumption at BS i. 

Power amplifier efficiency. 

Upper limit of the total energy consump¬ 
tion for BS i. 

Buying and selling energy prices; and 
functions thereof. 


C. Basic variables 


Vector channel from BS i to user k in slot t. 
iCf. Vector channel from all BSs to user k in slot t. 

\i\. Estimated channel for user k in slot t. 

Vector signal transmitted to user k in slot t. 
Information-bearing symbol of user k in slot t. 
yl, Received signal of user k in slot t. 

nj. Channel additive noise of user k in slot t. 

Bi Selection matrix forming the i-th BS’s transmit- 

beamforming vector. 

<5|, Channel estimation error of user k in slot t. 

y(, pi Proximal center and weight of the bundle method 
at iteration I. 

rji Predicted descent of the objective value for the 

bundle method at iteration 
gi^i Subgradient of G'(pi) at iteration i. 

r Energy buying-to-selling price ratio. 

E* Expected renewable generation. 

E* Realizations of random renewable generation. 

K, U Scaling variable and uniform random variable 
relating to E*. 


D. Decision variables 

Transmit beamforming vector from all BSs to user 
k in slot t. 
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Cl Stored battery energy at BS i at the beginning of 

slot t. 

Pi j Power delivered to or drawn from the batteries of 
BS i in slot t. 

Pg j Total energy consumption for BS i in slot t. 

Pi j Transmit power at BS i in slot t. 

El Energy harvested at BS i in slot t. 

Bi Vector collecting {El}f^i. 

Pi Auxiliary variable relating P^ ^ and 

Pi Vector collecting {P/}^i. 

X^. Matrix-lifting Beamformer for user k in slot t. 

tI Auxiliary variable introduced by the S-procedure. 

Z Matrix collecting all primal variables. 

A* Lagrange multiplier. 

A Matrix collecting all Lagrange multipliers. 


E. Functions 


SINRfc(.) 

SINR of user k. 

SINRfc(.) 

Worst-case SINR of user k. 

G(-,-) 

Worst transaction cost across entire horizon. 

G(.), G(.) 

Modified worst-case transaction cost. 

L(Z,A) 

Partial Lagrangian function. 

D{A) 

Dual function. 


I. Introduction 

To accommodate the explosive demand for wireless ser¬ 
vices, cellular systems are evolving into what are termed 
heterogeneous networks (HetNets) consisting of distributed 
macro/micro/pico base stations (BSs) to cover overlapping 
areas of different sizes m. Close proximity of many HetNet 
transmitters introduces severe inter-cell interference. Lor effi¬ 
cient interference management, coordinated multi-point pro¬ 
cessing (CoMP) has emerged as a promising technique for 
next-generation cellular networks such as LTE-Advanced El. 

To fully exploit the potential of CoMP at affordable over¬ 
head, coordinated beamforming and/or clustered BS coop¬ 
eration for downlink systems were investigated in ja-Ei. 
Multiple BSs cooperate to beamform, and each user’s data 
are only shared among a small number of BSs per cluster, 
thus greatly reducing the overall backhaul signalling cost. 

The rapid development of small cells in HetNets has also 
driven the need for energy-efficient transmissions. Due to the 
growing number of BSs, the electricity bill has become a 
major part of the operational expenditure of cellular operators, 
and cellular networks contribute a considerable portion of the 
global “carbon footprint” a. These economic and ecological 
concerns advocate a “green communication” solution, where 
BSs in cellular networks are powered by the electricity grid 
0-I1II1. However, the current grid infrastructure is on the 
verge of a major paradigm shift, migrating from the traditional 
electricity grid to the so-termed “smart grid” equipped with 
a large number of advanced smart meters and state-of-the- 
art communication and control links. To decrease greenhouse 
gas emissions, an important feature of future power systems 
is integration of renewable energy sources (RES). This leads 
to high penetration of distributed generators equipped with 
energy harvesting modules, which can crop energy from the 


environmental resources (e.g., solar and wind), and possibly 
trade the harvested energy with the main grid. In addition 
to distributed generation, distributed storage, and two-way 
energy trading associated with RES, demand-side management 
(DSM) including dynamic pricing and demand response, can 
further improve grid reliability and efficiency. Relying on 
pertinent tools, optimal energy management and scheduling 
with RES and/or DSM were proposed in ifT^ - lflSl . 

To take advantage of the aforementioned smart grid capa¬ 
bilities in next-generation cellular systems, only a few recent 
works have considered the smart-grid powered CoMP trans¬ 
missions ESl-lIISl. However, m only addressed dynamic 
pricing using a simplified smart grid level game, while ifTTl 
and ifTsl assumed that the energy amounts harvested from RES 
are precisely available a priori (e.g., through forecasting), and 
the harvested energy cannot be stored at the BSs. In addition, 
na assumed demand (or load) response based on different 
energy buying/selling prices across BSs without adapting 
power consumption to time-varying energy pricing. 

The present paper deals with optimal energy management 
and transmit-beamforming designs for the smart-grid powered, 
cluster-based CoMP downlink with the clustering carried by 
the HetNet’s central processor. Each BS has local RES and 
can perform two-way energy trading with the grid based on 
dynamic buying/selling prices. Different from ifTbll - lfTSll . we 
suppose that each BS has a local storage device, which can 
be charged to store the harvested (and even grid) energy and 
can be discharged to supply electricity if needed. To account 
for the stochastic and nondispatchable nature of both RES and 
wireless channels, we assume that the actual harvested energy 
amounts and the wireless channel states are unknown and time- 
varying, yet lie in some known uncertainty regions. Building 
on realistic models, we develop robust energy management and 
transmit-beamforming designs that minimize the worst-case 
energy cost subject to the worst-case user QoS guarantees for 
the CoMP downlink. Leveraging a novel formulation account¬ 
ing for the worst-case transaction cost with two-way energy 
trading, as well as the S-procedure in robust beamforming 
designs, we show how to (re-)formulate the said task as a 
convex problem. Strong duality of the latter is then utilized 
to develop a Lagrange dual based subgradient solver. It is 
shown that the resultant algorithm is guaranteed to find the 
desired optimal energy-management strategy and transmit- 
beamforming vectors, and could also facilitate distributed 
implementations among the BSs. 

The rest of the paper is organized as follows. The system 
models are described in Section II. The proposed approach to 
robust energy management and transmit-beamforming designs 
for the CoMP downlink is developed in Section III. Numerical 
results are provided in Section IV to demonstrate the proposed 
scheme. We conclude the paper in Section V. 

Notation: Boldface fonts denote vectors or matrices; 
is the K-hy-M dimensional complex space; (•)' denotes 
transpose, and conjugate transpose; tr(A) and rank(A) 
the trace and rank operators for matrix A, respectively; 
diag(ai,..., um) denotes a diagonal matrix with oi,..., om 
on the diagonal; | • | represents the magnitude of a complex 
scalar; A ^ 0 means that a square matrix A is positive semi- 
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Fig. 1. A two-BS CoMP system powered by smart grids, where BSs with local renewable energy harvesting and storage devices implement two-way energy 
trading with the main grid. 


definite. 


II. System Models 

Consider a smart-grid powered cluster-based CoMP down¬ 
link. A set I := of distributed BSs (e.g., 

macro/micro/pico BSs) provides service to a set K. := 
{1,... , A'} of mobile users; see Fig. [1] Each BS has M > 1 
transmit antennas and each user has a single receive antenna. 
Each BS is equipped with one or more energy harvesting de¬ 
vices (solar panels and/or wind turbines), and is also connected 
to the power grid with a two-way energy trading facility. 
Different from CSl-lIIl, each BS has an energy storage 
device (possibly consisting of several large-capacity batteries) 
so that it does not have to consume or sell all the harvested 
energy on the spot, but can save it for later use. 

Eor each CoMP cluster, there is a low-latency backhaul 
network connecting the set of BSs to a central controller 
0, which coordinates energy trading as well as cooperative 
communication. The central entity can collect both the com¬ 
munication data (transmit messages) from each of the BSs 
through the cellular backhaul links, as well as the energy 
information (energy buying/selling prices) from these BSs 
via smart meters installed at BSs, and the grid-deployed 
communication/control links connecting them. 

Assume slot-based transmissions from the BSs to the users. 
While the actual harvested energy amounts and wireless chan¬ 
nels cannot be accurately predicted, uncertainty regions for 
the wireless channels and renewable energy arrivals can be 
obtained, based on historical measurements and/or forecasting 
techniques. The slot duration is selected equal to the minimum 
time interval between the changes of the (channel or energy) 
uncertainty regions. Consider a finite scheduling horizon con¬ 
sisting of T slots, indexed by the set T := {1,... ,T}. Eor 


convenience, the slot duration is normalized to unity unless 
otherwise specified; thus, the terms “energy” and “power” will 
hereafter be used interchangeably throughout the paper. 


A. Downlink CoMP Transmission Model 

Per slot t, let G denote the vector channel from 

BS i to user k, Vi G I, Vfc € /C; and let := ,..., 

collect the channel vectors from all BSs to user k. With linear 
transmit beamforming performed across BSs, the vector signal 
transmitted to user k is 


ql = WfcSfc, vfc 


where denotes the information-bearing symbol, and G 
(j^M/xi denotes the beamforming vector across the BSs for 
user k. Eor convenience, is assumed to be a complex ran¬ 
dom variable with zero mean and unit variance. The received 
vector at user k is therefore 


yi = + 4 (1) 

l^k 

where is the desired signal for user k, ^^q? 

is the inter-user interference from the same cluster, and 
denotes additive noise, which may also include the downlink 
interference from other BSs outside this cluster. It is assumed 
that is a circularly symmetric complex Gaussian (CSCG) 
random variable with zero mean and variance cr^. 

The signal-to-interference-plus-noise-ratio (SINR) at user k 
can be expressed as 


SINRfe({w‘}) 
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The transmit power at each BS i clearly is given by 




(3) 


fee/C 


where the matrix g 

/ 


Bi := diag 


y (i-l)M M _ 


selects the corresponding rows out of to form the 

i-th BS’s transmit-beamforming vector. 

The channel state information is seldom precisely 
available a priori in practice. Relying on past channel mea¬ 
surements and/or reliable channel predictions, we adopt the 
following additive error model: = h^,where is the 

predicted channel. The uncertainty of this estimate is bounded 
by a spherical region mi: 


-f (5^ I ll^fell < e^l , Vfc,/ (4) 

where e]. > 0 specifies the radius of 7/^. This leads to the 
worst-case SINR per user k as [cf. @] 


SINRfe({w‘}) 


min - 




(5) 


To guarantee QoS per slot, it is required that 


^fc({w^})>7fc, Vfc (6) 


where jk denotes the target SINR value per user k. 


B. Energy Harvesting Model 

Let El denote the energy harvested during the last slot that 
is available at the beginning of slot t at each BS / € X; and let 
Gi := [E ]^..., EJ]'. Due to the unpredictable and intermittent 
nature of RES, is unknown a priori. In general, uncertain 
quantities can be modeled by postulating either an underlying 
probability distribution or an uncertainty region. Probability 
distributions (possibly mixed discrete/continuous) of the RES 
generation are seldom available in practice. Although (non- 
Iparametric approaches can be used to learn the distributions, 
the processes can be very complicated due to the spatio- 
temporal correlations incurred by various meteorological fac¬ 
tors. On the other hand, the approach of postulating an 
uncertainty region provides the decision maker range forecasts 
instead of point forecasts, which are essentially distribution- 
free. Suppose that lies in an uncertainty set £i, which can 
be obtained from historical measurements and/or fine forecast 
techniques. Erom the perspective of computational tractability, 
two practical paradigms for £i are considered here, 
i) The first model amounts to a polyhedral set mi; 

£!-.= le,\gi<El<Ei 

s 

T=\J 

t^7i,s s—1 




El 


< Ei:r, 



where {e\) denotes a lower (upper) bound on El', 
time horizon T is partitioned into consecutive but non¬ 
overlapping sub-horizons Ti^s, s = Each sub¬ 

horizon can consist of multiple time slots, and the total 
energy harvested by BS i over the sth sub-horizon is 
bounded by and 

ii) The second model relies on an ellipsoidal uncertainty set 
(see e.g. ||20l) 

£:r:={e, = + (8) 

where := [E^, ..., Ef]' denotes the nominal energy 
harvested at BS i, which can be the forecasted energy, or 
simply its expected value. Vector is the corresponding 
error in forecasting. The known matrix S 0 quantifies 
the shape of the ellipsoid £1, and hence determines the 
quality of the forecast. 

Remark 1: (Spatio-temporally correlated energy harvesting 
models). The aforementioned two practical models capture 
RES uncertainty across the scheduling (sub-)horizons per BS. 
The parameters required for constructing the sets {£’]’, } can 

be obtained offline via statistical learning techniques using 
historical data. In general, green energy harvested at different 
BSs can be spatially correlated if some BSs are geographically 
close. In this case, joint spatio-temporal uncertainty models 
can be postulated whenever the underlying correlations are 
known a priori', see details in mi- In general, a refined 
uncertainty model quantifying the actual harvested energy in 
a smaller region with a higher confidence level can be less 
conservative in the robust optimization formulation. However, 
the complexity of solving the resulting optimization problems 
directly depends on the choice of the uncertainty set. 

Note that the coherence time of RES arrivals can be much 
longer than that of wireless channels in practice ifTTll . IItTI . Yet, 
coherence times corresponding to the uncertainty regions of 
wireless channels can be much larger than those of the channel 
itself. In Section II-A, we implicitly assume that channel 
remains unchanged per slot. However, SINR^ in (|5]l can 
be easily redefined as the worst-case SINR over multiple 
channel coherence times with the same uncertainty region (and 
possibly different channel realizations). This way the issue of 
different time scales becomes less critical. In addition, our 
models in fact accommodate cases where uncertainty regions 
for RES arrivals remain unchanged over multiple time slots. 
The proposed approach presented next readily applies to obtain 
the robust ahead-of-time schedule in this setup. 


C. Energy Storage Model 

Let C° denote the initial energy, and C* the amount of 
stored energy in the batteries of the i-th BS at the beginning 
of time slot t. With bounding the capacity of batteries, 

it is clear that 0 < C* < Mi € X. Let denote the 

power delivered to or drawn from the batteries at slot t, which 
amounts to either charging (P^‘ > 0) or discharging (P^‘ < 
0). Hence, the stored energy obeys the dynamic equation 


Cl = C 


t-i 


Plr. 


t € 7”, % 


(9) 
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The amount of power (dis-)charged is also bounded by 


n, i s ^bi s ^b i 

< Pi 


b,i 


( 10 ) 


where < 0 and P™1^^ > 0 , while Wi € ( 0 , 1 ] is the 

battery efficiency at BS i. The constraint < Pl^ 

means that at most a fraction vJi of the stored energy is 
available for discharge. 


D. Energy Cost Model 

For the i-th BS per slot t, the total energy consumption 
Pi j includes the transmission-related power P^. j, and the rest 
that is due to other components such as air conditioning, data 
processor, and circuits, which can be generally modeled as a 
constant power. Pc ,, > 0 Ho), HU; namely. 


pli = Pc,^+pp/i = Po,^ + E 


fee/c 


where ^ > 0 denotes the power amplifier efficiency. For 
notational convenience, we absorb f into by redefining 
Bi := Bi/^; and further assume that Pj ^ is bounded by P™“. 

In addition to the harvested RES, the power grid can also 
supply the needed Pj ^ per BS i. With a two-way energy trad¬ 
ing facility, the BS can also sell its surplus energy to the grid 
at a fair price in order to reduce operational costs. Given the 
required energy P* j, the harvested energy P‘, and the battery 
charging energy Pl^, the shortage energy that needs to be 
purchased from the grid for BS i is clearly [Pj ^ — P* + P^ J ^; 
or, the surplus energy (when the harvested energy is abundant) 
that can be sold to the grid is [PP — P‘ -|- P^ J”, where 
[a]+ := max{a, 0}, and [a]“ := max{—a, 0}. Note that both 
the shortage energy and surplus energy are non-negative, and 
we have at most one of them be positive at any time t for 
BS i. 

Suppose that the energy can be purchased from the grid 
at price a‘, while the energy is sold to the grid at price /?‘ 
per slot t. Notwithstanding, we shall always set > /3‘, V, to 
avoid meaningless buy-and-sell activities of the BSs for profit. 
Assuming that the prices a‘ and /3‘ are known a priori, the 
worst-case transaction cost for BS i for the whole scheduling 
horizon is therefore given by 

T 

G{{PP}, {Pi,,}) := max ^ 

‘ i=l 

-/3‘[PP-P* + P,P]-)- (11) 


III. Energy Management for CoMP Beamforming 

Based on the models of Section [III we consider here robust 
energy management for transmit beamforming in a CoMP 
cluster. Knowing only the uncertainty regions of the wireless 
channels and renewable energy arrivals, the central controller 
per cluster performs an (e.g. hour-) ahead-of-time schedule 
to optimize cooperative transmit beamforming vectors {w^} 
and battery charging energy {P^ j}, in order to minimize the 
worst-case total cost G{{Pli}, {Pl,i}), while satisfying 


the QoS guarantees SINRfe({w^}) > 7 ^, Vfc, over the schedul¬ 
ing horizon T- Eor convenience, we introduce the auxiliary 
variables Pi := PP -f Pl^^, and formulate the problem as 

minimize G({P/}) (12a) 

suject to; 

p ! = Pc,, + E + p,p, Vi e I, f G r (12b) 

keK 

0 < Pc,, + E < Plir, (12c) 


k&K 

Cl = cl-^ + Pi,, yiGi,t€r (i 2 d) 

0 < G- < Vi G X, f G r (12e) 

< pi,, < p^r, vi g x, f g r (Uf) 

- w,cl-^ < Pl„ Vi G X, f G r (12g) 

^fc({w^})> 7 fc, ykG]C,tGT. ( 12 h) 


Consider for simplicity that (fT2l i is feasible. The problem can 
in fact become infeasible if the SINR thresholds 7 ^ are high 
and the wireless channel qualities are not good enough. In 
this case, the subsequent problems (I22al i are always infeasible. 
Recall that the proposed scheme is proposed to determine 
the ahead-of-time beamformer and energy schedule (offline). 
Such an infeasibility, once detected, can naturally lead to an 
admission control policy, i.e., to a criterion for dropping users 
or SINR requirements that render the problem infeasible ll22l . 

Solving (fT^ can provide a robust solution for the smart-grid 
powered CoMP downlink with worst-case performance guar¬ 
antees. It is worth mentioning here that thanks to the worst- 
case cost G({P/}), randomness introduced due to the wireless 
fading propagation and also due to the RES uncertainty 
has been eliminated; thus (fT 2 T i contains only deterministic 
variables. Because of (I12bb . (I12cl i. and (I12hb . the problem is 
nonconvex, which motivates the reformulation pursued next. 


A. Convex Reformulation 

Eirst, let us consider convexity issues of the objective 
function G({P/}). Define := (a* — I3'‘')j2 and (jf := 
(a‘ -|- /?‘)/2, and then rewrite: 

T 

G({i^‘}) = max E {^*\Pl - El\ + 4 >\pI - P‘)) • 

‘ t=i 

Since > /3‘ > 0, we have cjf > > 0. It is then clear 

that I P-‘ — P‘ I -I- (Pi—El) is a convex function of P-‘ for 
any given Ef As a pointwise maximization of these convex 
functions, G({P/}) is also a convex function of {P/} (even 
when the set £i is non-convex) ll2Tl . 

Except for (I12bb . (I12cb . and (I12hb . all other constraints 
are convex. We next rely on the popular semidefinite pro¬ 
gram (SDP) relaxation technkju^ to convexify (I12hb . By 
the definitions of Hi and SINRfe({w^,}), the constraint 
SINRfe({w^}) > 7 fe can be rewritten as; 

Fki^l) > 0 for all 6} such that dl < (e^)^ (13) 
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where 


F,idl):= {hl + diri 


li 


l^k 


Using ( foi l and upon applying the well-known S-procedure 
in robust optimization ll24l . the original problem (fT2l i can be 
reformulated as an SDP with rank constraints. 

(S-procedure) Let A,B be n x n Hermitian matrices, c S 
C" and d €M. for which the interior-point condition holds; that 
is, there exists an x such that x^Bx < 1; and the following 
are equivalent: 

(i) x^Ax -f- c^x -f x^c + d > 0 for all x S C" such that 
x^Bx < 1; 

(ii) There exists a z/ > 0 such that 


A + i/B 


C 

d — V 


>- 0 . 


For := S it clearly holds that X*^ ^ 

0 and rank(X^) = 1. Using the S-procedure, (fOl) can be 
transformed to 


r 


t 

k 


fn+rll 


■\rt ut 
* k’^^k 


^ k'-'-k ^k 



> 0 


(14) 


where > 0 and 

Y‘:=^X‘-^X*. (15) 

i^k 


Introducing auxiliary variables and dropping the rank 
constraints rank(X^) = 1, \/k,t, we can then relax (fT^ to: 


(16a) 

(16b) 
(16c) 
(16d) 

In addition to the linear constraints ( I12db -( [T2^ , the quadratic 
power constraints (I12bb and ( fna have now become linear, 
and SINR constraints (I12hb become a set of convex SDP 
constraints in (HSll. Since the objective function is convex, 
problem (fTbl l is a convex program that can be tackled by a 
centralized solver; e.g., using the projected subgradient descent 
approach. However, the feasible set (I16bb - (ll6db is the inter¬ 
section of a semi-definite cone and a polytope, for which the 
iterative projection is complicated. To reduce computational 
complexity and enhance resilience to failures, we next develop 
an efficient algorithm to solve (fTbl) in a distributed fashion 
coordinated by different agents. 



B. Lagrange Dual Approach 

Since dlS} is convex, a Lagrange dual approach can be 
developed to efficiently find its solution. Let A*, denote 
the Lagrange multipliers associated with the constraints (I16bb . 


With the convenient notation Z := {X^, r^., C‘, P/} and 

A := {A*}, the partial Lagrangian function of (fTbl l is 

P(Z,A) := 

T 

E+E + E MB,x^)+P‘, - pi ). 

iei t=i keK 

(17) 


The Lagrange dual function is then given by 
D{A) ■= mm P(Z, A) 
s. t. ([Hdli - ([T^, (flScli - (fTM 
and the dual problem of (fTbl l is: 


( 18 ) 


max D{A). (19) 

A 

Subgradient iterations: Let j denote the iteration index. To 
obtain the optimal solution A* to the dual problem (fT^ . we 
resort to the dual subgradient ascent method, which amounts 
to the following update 


>Hij + 1) = >hU) + (20) 

where /i(j) > 0 is an appropriate stepsize. The subgradient 
§( 7 ) ■= bA*(j) can then be expressed as 

=i"c,* + Etr(B*XUj))+P(,(j)-i^‘(j), Vz,f (21) 

keK 

where X^(j), Pi ,^{ 3 ), and P/(j) are given by 


{Xfc(i)}Li Garg min El^iO') E 

s. t. (I16cb — (I16db . VA: (22a) 

T 

{A*z(7)}Li e arg min Yl[Xl{j)Pli] 


s. t. (I12db — ( |12g| i, Vf 


(22b) 


{^*(7)}Li e arg mm [G{{Pl}) - E • (22c) 


Solving the subproblems: Subproblems (I22ab are standard 
SDPs per t G T; hence, {X^(j)}^j^ for all t can be efficiently 
solved by general interior-point methods ll23l . 

The subproblems (I22bl i are simple linear programs (LPs) 
over {P‘,, per i G X; hence, {Pl,U)}I=i^ Vt, can be 

obtained by existing efficient LP solvers. 

Due to the convexity of G({P,‘}), the subproblems (I22cb are 
convex per i G X. Yet, because of the non-differentiability of 
G({P,‘}) due to the absolute value operator and the maximiza¬ 
tion over e, G £i, the problem is challenging to be handled 
by existing general solvers. For this reason, we resort to the 
proximal bundle method to obtain {P,‘(j)}^]^. Upon defining 
G(p0 := G(pO-ELi X\{j)Pl, where pi := [P^, ... 
the partial subgradient of G(pi) with respect to P,‘ can be 
obtained as 


t)G(p.) \a\-\\{3), ifP/PPf 

dPl \/3|-AKj), ^PI<eT 




















7 


where e* := , Ef*]' for the given is obtained as 

T 

e* G arg max ^ \P* - El\ + {P* -El)). (24) 

‘ t=i 

It can be readily checked that the objective function in (l24l i is 
convex in under the condition a* > [3*, € T. Hence, the 

globally optimal solution is attainable at the extreme points 
of Ei- Leveraging the special structure of Ei, we utilize an 
efficient vertex enumerating algorithms to obtain e* directly, 
as detailed next. 


C. Proximal Bundle Method 

Given the partial subgradient in (l27t . nonsmooth convex 
optimization algorithms can be employed to solve the subprob¬ 
lem (I22cb . A state-of-the-art bundle method ||2^ Ch. 6 ], Il25l 
will be developed here with guaranteed convergence to the 
optimal see also ifTTl . 

Similar to cutting plane methods, the idea of bundle method 
is to approximate the epigraph of an objective by the inter¬ 
section of a number of halfspaces. The latter are generated 
through the supporting hyperplanes, referred to as cuts, by us¬ 
ing the subgradients. Specifically, letting i denote the iteration 
index of the bundle method, the iterate Pi/+i is obtained by 
minimizing a polyhedral (piecewise linear) approximation of 
G{pi) with a quadratic proximal regularizer 

:= argminjc^jp,) -f ^||p, - | (25) 

PiSR'^ t ^ J 

where := max{G(p*,o) + go(p*-P*,o), ■ ■ •, G'(p*.^) + 

“ Pi,()}’ is the subgradient of G'(pi) evaluated at 
the point p = [cf. (I23ll1; and the proximity weight 
controls the stability of iterates. 

Different from the proximal cutting plane method (CPM), 
the bundle method updates its proximal center yg according 
to a descent query 


yt+i = 


Pi,t+i, if G(y^) - G(pi,^+i) > % 
yt, if G{yi) - G(p*,f+i) < % 


(26) 


where 0 G (0,1) is a pre-selected constant, and 77 ^ := G{y() — 
^Gf(p£+i) + ^||p^_|_i — y^lP^ is the predicted descent of the 
objective in (IZSl l. Essentially, if the actual descent amount 
G{yt) — G(p^+i) is no less than a 0 fraction of the predicted 
counterpart then the iterate takes a “serious” step updating 
its proximal center y^+i to the latest point p^+i; otherwise it 
is just a “null” step with the center unchanged. The intelligent 
query (l26l l enables the bundle method to find “good” proximal 
centers along the iterates, and hence it converges faster than the 
proximal CPM. In addition, depending on whether a serious 
or a null step is taken, the proximity weight can be updated 
accordingly to further accelerate convergence 12 ^ : that is, 

^ f max(pf/10,pmin), if G(y£) - G(pi,^+i) > % 

1 min(10p^,pmax), if G(y^) - G(pi,^+i) < %. 

The algorithm terminates if y^ = Pi,^+i, while finite termina¬ 
tion is achievable if both the objective and the feasible set are 
polyhedral ||23] Sec. 6.5.3]. 


Now, to complete the proximal bundle method for solving 
(I22cb . we only need solve problem dZST l. Using an auxiliary 
variable u, dZSl) can be re-written as 

min u + ^\\p, - yiW^ (27a) 

Pi,u Z 

s.t. G(pi,„)-|-g'_„(pi - < u, n = 0,1,... ,f. (27b) 

Introducing multipliers ^ G and letting 1 denote the 

all-ones vector, the Lagrangian function is given as 

£{u,pi,^) = (1 - l'|)u+ y||p, - y^f 
i 

+ X! + Sl.uiPi - P»,n))- (28) 

n—0 

The optimality condition Vp.C{u,Pi,^) = 0 yields 


Pi — y( ^ ' ^nSi,n ■ 

Pi n 

n—0 

Substituting < |29l t into < |28] t, the dual of < |27] t is 


max- 

4 


n—0 

s.t. 1 ^ 0 , I'l = 1 . 


^in{G{p,^n) + ^i,n{yi 


n—0 


(29) 


Pi,n)) 

(30a) 

(30b) 


It can be readily seen that (l30l l is essentially a QP over the 
probability simplex in the dual space, which can be solved 
very efficiently as in e.g., M- 


D. Optimality and Distributed Implementation 

When a constant stepsize p,{j) = p is adopted, the 
subgradient iterations (|20ll are guaranteed to converge to a 
neighborhood of the optimal A* for the dual problem (fT^ 
from any initial A(0). The size of the neighborhood is 
proportional to the stepsize p. If we adopt a sequence of non- 
summable diminishing stepsizes satisfying limj_>oo pU) = 0 
and ~ iterations (l20b will asymptoti¬ 

cally converge to the exact A* as j ^ 00 ||23]| . 

The objective function (I16ab is not strictly convex because 
it does not involve all optimization variables. Hence, when 
it comes to primal convergence, extra care is necessary ll28ll . 
II 29 I . Specifically, the optimal primal can be attained either by 
adding a strictly convex regularization term, or, by utilizing 
the augmented Lagrangian. Here, we will simply implement 
the Cesdro averaging method ll30l to obtain the optimal power 
schedules. With p'^^ := m(j)> the running average is 

- m 

(31) 

Msum 

which can be efficiently computed in a recursive way 

jrn . (32) 

,,m ^ ' i,m 

A^sum A^sum 

Note that if a constant stepsize p{j) = p is adopted, OTI) and 
(I 32 I) boil down to the ordinary running average 

1 ™ 1 — 1 

Z"* = — y Z{j) = —Z(m) -b —— -Z^-^. 
m m m 










TABLE I 

Generating capacities, battery initial energy and capacity, 

CHARGING LIMITS AND EFEICIENCY. 


Unit 

pmin 

pmax 

^Gi 


^max 

pmin 

pmax 

^bi 

TXJ'l 

1 

0 

50 

5 

30 

-10 

10 

0.95 

2 

0 

45 

5 

30 

-10 

10 

0.95 

3 

0 

45 

5 

30 

-10 

10 

0.95 

4 

0 

45 

5 

30 

-10 

10 

0.95 

5 

0 

50 

5 

30 

-10 

10 

0.95 

6 

0 

45 

5 

30 

-10 

10 

0.95 


TABLE II 

Limits of forecasted wind power and energy purchase prices 


Slot 

1 

2 

3 

4 

5 

6 

7 

8 

Si 

2.47 

2.27 

2.18 

1.97 

2.28 

2.66 

3.10 

3.38 

si 

2.57 

1.88 

2.16 

1.56 

1.95 

3.07 

3.44 

3.11 

si 

2.32 

2.43 

1.27 

1.39 

2.14 

1.98 

2.68 

4.04 

si 

2.04 

1.92 

2.33 

2.07 

2.13 

2.36 

3.13 

4.16 

s| 

2.11 

1.19 

2.26 

2.19 

1.55 

2.71 

3.37 

2.45 

si 

2.01 

2.29 

2.20 

0.98 

2.43 

3.22 

2.74 

3.93 


0.402 

0.44 

0.724 

1.32 

1.166 

0.798 

0.506 

0.468 


If for the relaxed problem (HSll, the obtained solution 
satisfies the condition rank(X^*) = 1, then it clearly 

yields the optimal beamforming vectors w^* as the (scaled) 
eigenvector with respect to the only positive eigenvalue of 
X^* for the original problem (fT2l i. Fortunately, it was shown 
in lISTl Them. 1] that the S-procedure based SDP (fThl l always 
returns a rank-one optimal solution X^*, when the 

uncertainty bounds are sufficiently small. If Ck is large, 
the existence of rank-one optimal solutions of (fThl l cannot 
be provably guaranteed. In this case, a randomized rounding 
strategy needs to be adopted to obtain vectors from 
X^* to nicely approximate the solution of the original problem 
d. Even though no proof is available to ensure a rank-one 
solution when Cfc is large, it has been extensively observed in 
simulations that the SDP relaxation always returns a rank-one 
optimal solution OTIl . This confirms the view that the optimal 
beamforming vectors for the original problem (fTSli will be 
obtained by our approach with high probability. 

The subgradient iterations can be run in a distributed 
fashion. Specifically, the central controller can maintain the 
Lagrange multipliers A(j) and broadcast them to all BSs via 
backhaul links. Given the current A(j), the central controller 
solves the subproblems (122al) to obtain the beamforming 
vectors for all BSs. On the other hand, each BS solves its 
own subproblems (122bl i- (122cl) . which are decoupled across 
BSs. The BSs send back to the central controller {P* 

{Pb {PiU)}J=i’ which are in turn used to 

update A(j + 1) through the subgradient iterations (l20l) . 

IV. Numerical Tests 

In this section, simulated tests are presented to verify 
the performance of the proposed approach. The Matlab- 
based modeling package CVX 2.1 13^ along with the 
solvers MOSEK 7.0 041 and Sedumi 1.02 051 are used 
to specify and solve the resulting optimization problems. All 
numerical tests are implemented on a computer workstation 
with Intel cores 3.40 GHz and 32 GB RAM. 



Iteration 



Fig. 2. Convergence of the objective value and the subgradient norm of 
Lagrange multipliers (M = 2, K = 10, r = 1). 



Iteration 


Fig. 3. Convergence of the bundle method solving the subproblem (Hi 
{M = 2, K = 10, r = 1, j = 2). 


The scheduling horizons of the considered CoMP network 
is T = 8. Two configurations are tested: (Cl) 1 = 2 BSs 
and K = 10 end users (small size); and (C2) 1 = 6 BSs 
and K = 30 end users (large size). All wireless channels 
are assumed to be flat Rayleigh fading, and normalized to 
unit power. The noises are modeled as circularly-symmetric 
Gaussian random vectors. Without loss of generality, the ef¬ 
fects of path loss, shadowing, and Doppler fading are ignored. 
Parameters including the limits of Pc^, Ci, Pbi and the 
discharging efficiency Wi are listed in Table |T] A polyhedral 
uncertainty set (l7]i with a single sub-horizon (no time partition) 
is considered for the RES. In Table |II] the energy purchase 
price a* is given across the entire time horizon. The selling 
price is set as /3‘ = ra* with r € [0, 1]. In addition, the 
lower limits are listed therein, which were rescaled 

from real data we obtained from the Midcontinent Independent 
System Operator (MISO) ll36l . The upper limits were set to 
Pi = while the total horizon bounds are = 

El- 

First, convergence of the objective value (116ab . and the 
^ 2 -norm of the subgradient of the running-average Lagrange 
multiplier (EB is verified in Eig. ID It can be seen that 
both metrics converge within a few hundred iterations, which 
confirms the validity of the dual decomposition approach 
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Fig. 4. Optimal power schedules of P* (I = 2, M = 2, K = 10, r = 1). Optimal power schedule of {I — 2, M — 2, K — 10, r — 1). 




Fig. 5. Optimal power schedules of P/ (I = 6, M = 2, K = 20, r = 1). Fig. 7. Optimal power schedule of P^. {I = 0, M = 2, K = 20, r = 1). 


along with the subgradient ascent algorithm. With the Ces^o 
averages, convergence of the dual and primal variables was 
also confirmed, but it is omitted due to limited space. 

Figure [3 depicts the effectiveness of the proposed bundle 
method minimizing the convex nonsmooth objective (I22cl l. 
Clearly, incorporating the scheme of dynamically changing the 
proximity weight pi, the bundle algorithm converges very fast; 
typically, within 10 iterations. 

The optimal power schedules of P/ are depicted in Figs. |4] 
and 12 For both configurations, the stairstep curves show that 
the lowest levels of P/ occur from slot 4 to 6 . This is because 
the energy selling and purchase prices are relatively high 
during these horizons [cf. Sec. m, which drives the BSs’s 
power consumption low in order to minimize the transaction 
cost. 

Figs. m and | 2 ] show the optimal battery (dis-)charging 
amount Pg , while Figs. 0 and |9] depict the state of charge 
C*. As a component part of P*, P|j exhibits a similar trend 
in response to the price fluctuation; that is, there is a relative 
large amount of battery discharging (Pg. < 0) when the 
corresponding price is high. Note that both Pg and C* never 
exceed their lower and upper limits [cf. Table HI . 

Robustness of the worst-case design to the uncertainty of 
channel estimates [cf. @] is confirmed in Figs. [TOl and [TT] 
The red solid line indicates the SINR threshold 7 ^ = 0.1 


that is set common to all users for simplicity. The non-robust 
scheme simply considers the estimated channel as the 
actual one, and plugs it into the worst-case SINR control 
design (| 6 ]l. This constraint can be relaxed to a linear matrix 
inequality: ~ which is a relaxed version 

of our proposed counterpart in (O. For both the robust 
and non-robust approaches, each transmit beamformer 
is obtained as the principal eigenvector associated with the 
largest eigenvalue of the resulting X^. The CDF of the actual 
SINR is obtained by evaluating (|2]l with 5, 000 independent 
and identically distributed (i.i.d.) channel realizations. The 
channel perturbations are first generated as complex 

Gaussian, and then rescaled to the boundary of the spherical 
region [cf. ([ 111 ]. It can be seen that 20% of the realizations 
of the non-robust scheme does not satisfy the SINR constraint, 
while only about 2% for the proposed approach. Note that the 
SDP relaxation is not always exact, which results in violating 
the SINR threshold for a few channel realizations. 

Finally, CDFs of the transaction cost are depicted in Fig. [12] 
The proposed robust approach is compared with a heuris¬ 
tic scheme that assumes the expected renewable generation 
E* = + Ej^) is available and prefixed for problem ( [T2 ). 

The CDF curves were plotted by evaluating the transaction 
cost ( IIbal l with 10 ® realizations of t and the obtained 

optimal solutions {.P/ji,*. The RES realizations were gener- 
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Time slot 


Fig. 8. Optimal power schedule for Cj (I = 2, M = 2, K = 10, r = 1). 



Time slot 


Fig. 9. Optimal power schedule for Cj {I = 0, M = 2, K = 20, r = 1). 

ated as {El}i,t = Md + i^U{e\ — ^*), where is a uniform 
random variable on [0,1]. Three cases with k = 0.01,0.1, 0.5 
were tested. Clearly, transaction costs of both the proposed 
and the heuristic methods decrease with the increase of k. 
Since a larger value of k implies more harvested renewables 
energy yielding a reduced transaction cost. Note that negative 
transaction costs means net profits are obtained by selling 
surplus renewables back to the smart grid. Interestingly, for 
all cases, the proposed approach always outperforms the 
heuristic scheme with less transaction costs. This is because 
the proposed schedules of {Pi i,Pli}i,t are robust to the 
worst-case renewable generation {El}i^f However, in practice 
more RES is often available than the worst case. Hence, the 
proposed method has a larger profit-making capability than the 
heuristic alternative. 

V. Conclusions 

Robust ahead-of-time energy management and transmit- 
beamforming designs were developed to minimize the worst- 
case energy costs subject to the worst-case user QoS guar¬ 
antees for the CoMP downlink, which is powered by a grid 
with smart-meter based dynamic pricing and RES available 
at the BSs. Building on innovative models, the task was 
formulated as a convex program. Relying on a Lagrange- 
dual approach, efficient algorithms were introduced to obtain 



SINR 


Fig. 10. SINR cumulative distribution function (CDF) {I = 2, M = 2, 
K = 10,r = 1). 



Fig. 11. SINR cumulative distribution function (CDF) (/ = 6, M = 2, 
K = 20,r = 1). 

optimal solutions. The proposed scheme provides the offline 
ahead-of-time beamformer and energy schedule over a finite 
time horizon. Development of online management schemes 
will be an interesting direction to pursue in our future work. 

Supported by major programs such as the US Energy 
Independence and Security Act and European SmartGrids 
Technology Platform, the smart grid industry has seen fast 
growth. It is expected that all future communication systems 
will be powered by a smart grid. As a result, integration 
of smart-grid technologies into system designs will hold the 
key to fully leverage energy-efficient communications in next- 
generation HetNets. The proposed models and approaches 
can pave a way to further advancing fundamental research 
on smart-grid powered HetNet transmissions, which will be 
pursued in future works. 
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